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Abstract 

We introduce a new set of algorithms to compute Jacobi matrices associated with mea- 
sures generated by infinite systems of iterated functions. We demonstrate their relevance in 
the study of theoretical problems, such as the continuity of these measures and the logarith- 
mic capacity of their support. Since our approach is based on a reversible transformation 
between pairs of Jacobi matrices, we also discuss its application to an inverse / approxima- 
tion problem. Numerical experiments show that the proposed algorithms are stable and can 
reliably compute Jacobi matrices of large order. 
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1 Introduction 

Let /X be a non-negative probability measure with compact support in R and let J^^ be the 
associated Jacobi matrix: 

/ ao bi \ 
bi ai 62 

62 a2 63 ■ 



V •■ / 

This infinite symmetric tridiagonal matrix encodes the three-terms recurrence relation of the 
orthonormal polynomials {pn{fJ.; s)}neN of /i: 

sp„(^;s) = 5„+ip„+i(/x; s) -|-a„p„(^;s) -I- 6„p„_i(/i; s), (2) 

with 60 = and initialized by s) — 0, Po(m; s) — 1. Recall that the integral / Pn(/i; s)pm{p] s)dfi{s) 

is equal to one when n = m and is null in all other cases. 

According to Gautschi [23], the computation of the Jacobi matrix associated with a given 
measure is a fundamental problem of numerical analysis. In fact, when the moment problem is 
determined [2] (this is the case under the compactness hypothesis above), the Jacobi matrix 
uniquely identifies the measure fi. Its computation is necessary for the evaluation of orthogonal 
sums via Clenshaw's algorithm |9j , but arguably its most important role comes in differentiation 
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and integration using orthogonal polynomials |14j and in Gaussian quadratures [241 125[ I36j , via 
the linear algebra technique of Golub and Welsch [211 130] • 

Jacobi matrices play a major role also in mathematical physics, since they can be seen as 
discrete Schrodinger operators acting in ^2(Z+), properly when 6„ = 1 for all n and in an extended 
sense in the general case. The investigation of the links between the properties of the sequences 
of coefficients a„ and 6„ and those of the "spectral" measure fj, has involved both students 
of orthogonal polynomials and numerical analysis and of mathematical physics and quantum 
dynamics, see e.g. [51], [10] and references therein. In the latter setting the Fourier transform 
of jj, and of its orthogonal polynomials take on the dynamical meaning of probability amplitudes 
of the quantum motion [461 1471 148j . Remark that they can be numerically computed in a very 
efficient way directly from the Jacobi matrix [45j . 

The asymptotic properties of the Fourier transform of a measure have been the focus of 
intense investigation also in harmonic analysis [SS] HT] [551 ISS] : these studies have highlighted 
the distinctive properties of singular continuous measures, which appear rarely on the stage of 
numerical analysis, but that, to the contrary, are principal in this paper. 

In fact, we consider herein measures generated by Systems of Iterated Functions (I.F.S) [34] l3]. 
that are frequently, although not always, singular with respect to the Lebesgue measure. I.F.S. 
are a highly versatile tool, that has been employed for the approximation of natural objects as well 
as for image compression [T3l|5], for wavelet construction [15] and for numerical integration [7]. In 
this forcefully minimal list of applications of I.F.S. an item of theoretical relevance must be added: 
Jacobi matrices of I.F.S. measures have been used as abstract models of aperiodic Schrodinger 
operators [32l |42l |43l , whose fine spectral properties are crucial for the phenomenon of wave 
propagation, in a study which lies at the intersection of the three disciplines just mentioned: 
numerical analysis, mathematical physics and harmonic analysis. 

Computing the Jacobi matrix associated with an I.F.S. measure is therefore a fundamental 
task, but a challenging one, the more so because the usual algorithms based on modified moments 
[55] are ill conditioned, for reasons explained in [HI [41] (see also |6]). This goal has been 
achieved with ad hoc techniques [HJ [191 Hll Hi] for I.F.S. with finitely many maps. In parallel, 
the computation of the Jacobi matrix for refinable functionals [501 HZ] (a particular case of 
I.F.S.) provided results [38l [39] that can compared to those for I.F.S. Yet, not all of the above 
algorithms are capable of producing Jacobi matrices of large orders, a crucial requirement for the 
investigation of the theoretical questions mentioned before. 

In this paper, we contribute a new set of algorithms to this family, that are numerically stable 
to large orders and that are also capable of handling the case of I.F.S. with uncountably many 
affine maps. The proposed algorithms are based on the existence of a one-to-one transformation 
between the Jacobi matrix of /i and that of an auxiliary measure, cr, that encodes the parameters 
of affine, homogeneous I.F.S., that we will define momentarily. 

Reversibility of such transformation also enables us to solve an inverse / approximation 
problem — that of computing the set of I.F.S. maps generating a given measure /i. Original results 
[7] indicate that I.F.S. quadratures, derived from the solution of the inverse problem, might offer 
advantages with respect to conventional ones — yet, as was to be expected, numerical stability is 
an important issue also in this inverse problem. Despite previous efforts [l] [5l [71 [33l [1] [21] the 
first algorithm really achieving large orders in a stable way has been announced in [49]: since it 
completes the theory developed herein and it can serve to prove experimentally the stability of 
the forward algorithms, it will be briefly discussed in the final section of this paper. 

This paper is organized as follows: in the next section we review the formalism of homogeneous 
iterated function systems and we generalize it to allow for uncountably many maps. This leads to 
the definition of a convolution-like operator on measures, that is studied in Sect. [S] particularly 
with respect to its action on orthogonal polynomials and Jacobi matrices. This section is the 
core of the paper, since it contains almost all technical lemmas. In Sect. HI we use these lemmas 
to build numerical algorithms for the computation of the I.F.S. convolution and of the Jacobi 
matrix associated with an I.F.S., that are experimentally examined with respect to stability and 
performance in a number of significant cases. In Sect. Othese properties and the preceding theory 
allow us to derive rigorous results on, as well as numerical estimates of, significant analytical 
properties of I.F.S. measures such as their continuity and the capacity of their support. In Sect. 
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[6]we describe a different technique for computing I.F.S. convolutions and Jacobi matrices, based 
on discrete measures and on inverse Gaussian methods. To complete the paper, in Sect. [7] we 
discuss the inverse algorithm, that yields the I.F.S. approximations of a target measure. The 
conclusions, Sect. [H briefly sum up the work and discuss further developments. 

2 Homogeneous AfRne Iterated Function Systems 

Systems of iterated functions |34j[3l[5] are finite collections of real maps : R — )• R, i = 1, . . . , M, 
for which there exists a set A, called the attractor of the I.F.S., that solves the equation A = 
Ui=i M 4>i{^)- Existence and uniqueness of A can be ensured under general circumstances. 
In this paper, we adopt a specific choice of the maps (j) that has the advantage of leading to a 
structured algebraic problem: that of homogeneous, affine transformations: 

Us) = Sis - + /3, i = l,...,M, (3) 

where (5 is a real constant between zero and one, and (3i are real constants, that geometrically 
correspond to the fixed points of the maps. By associating a positive weight, tt^ > 0, tt^ = 1, 
to each map one can define a measure fi supported on A via a procedure that we generalize in eq. 
([5]) below. The measure /i is specified by the value of S and of the pairs (/3i, tt^), for i — 1, . . . , M. 
This completes the definition of a "classical" homogeneous I.F.S. 

Rather than restricting the cardinality of maps, M, to be finite (or countable, as done by 
Maulin and Urbansky |51| ) we now allow the index set to be a continuum. This can be done in 
a variety of ways |52j . Our approach is to follow \16\ and to observe that a finite, homogeneous 
I.F.S. is fully described by the choice of S and of the discrete measure 

M 

where Dr^ is a unit mass atomic (Dirac) measure at the point x. We now let any positive 
probability measure a to be the distribution of affine constants: we only assume that the support 
of (T is contained in a finite interval, which, without loss of generality, we may take to be [—1,1]. 

Definition 1 Let a be a positive Borel probability measure on R whose support is contained in 
[—1, 1], let S be a real number in [0, 1) and let S := 1 — 6 . Let the real number /3 parameterize the 
LF.S. maps 4'{I3, ■) as 4>{p, s) := Ss + 5/3. The invariant LF.S. measure associated with the affine 
homogeneous I.F.S. (6, a) is the unique probability measure /i that satisfies 

J f{s) dfi{s) ^ J da{p) j dfi{s) f{Hl3,s)), (5) 

for any continuous function f . 

Remark 1 Lf 6 — 0, then ii ^ a. Ln fact, in this case 4>{(3, s) = /3 for all s. 

Consistency of Def. [1] i.e. existence and uniqueness of /i are easily proven to hold: 

Proposition 1 For any 6 £ [0, 1) eq. defines an invertible transformation from the space 
A4([— 1, 1]) of probability measures a on [—1, 1] to the subset o/A^([— 1, 1]) composed of invariant 
measures ji of homogeneous LF.S. with contraction ratio S. The support of ii contains the support 
of a and the convex hulls of these two sets coincide. 

Proof. This result is already contained in [16), although not entirely, because of a different 
parametrization of I.F.S. maps. Let us start from the second statement. Recall that the attractor 
A of an I.F.S., which is the support of the invariant measure /i, can also be characterized as the 
closure of the set of fixed points of the composition of an arbitrary number of maps 4>{P, ■). It 
is then apparent that the support of cr is a subset of that of fi, and that the convex hull of the 
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support of fi is the interval between the infimum and the supremum of the set of fixed points l3, 
that is, the convex hull of the support of a. 

Next, since the supports of a and /x are contained in a finite interval, their moment problem is 
determined [2] and their infinite sets of moments uniquely identify the corresponding measures. 
The two sets of moments can be put in one-to-one relation: putting /(s) = in eq. ([5]), for any 
j, one gets a triangular set of relations by which it is possible to derive all moments of /i from 
those of a, and viceversa. I 

Remark 2 The proof above does not offer a numerically viable mean of computing the described 
correspondence between n and a, at least as far as Jacobi matrices are concerned, as noted in the 
Introduction. 

3 A convolution-like operator and orthogonal polynomials 

We now want to define a mapping $ in the space A^([— 1, 1]) of probability measures on [—1,1], 
according to which the invariant measure of an I.F.S. is the fixed point of such transformation: 
II = Since $ turns out to be contractive in a suitable metric, /i can be found as the limit of 

the sequence of measures /i„ :— $"/io, where /iq is any initial probability measure. To achieve this 
goal, we first describe a convolution-like operator of measures induced by the I.F.S. construction 
and derive its basic analytical properties (next subsection). We then describe the scaling relations 
of orthogonal polynomials with respect to afRne transformations (second subsection). Thanks to 
these relations, the action of the operator <i> can be finally transferred on the Jacobi matrices of 
measures in A^([— 1, 1]) (third subsection). 

3.1 I.F.S. convolution of measures 

We define a convolution like-operator ^s{'', ■) acting on A^^([— 1, 1]) as follows. 

Definition 2 Let a, rj two measures in A^([— 1, 1]) and let 5 G [0, 1). The measure fj e A^([— 1, 1]), 
called the IPS convolution of a and rj, fj := <^s[a'.Tii), is defined via the equation 

J f{s) dfj{s) = J da{p) j d7is) fm,s)), (6) 

holding for any continuous function f . 

Remark 3 The symmetric role of a and rj in eq. (0j shows that they can be interchanged in the 
action of ^s{-; •), provided one also exchanges 6 and S: ^^(cr;?;) = ^g(ri;a). 

Suppose now that one keeps a fixed in the above construction. Then, fj := $5 (cr; 77) is a 
function of 77 alone: this defines the mapping $5(0'; •) in A4([— 1, 1]). According to well established 
theory, see e.g. \W\ I52|. one can easily prove 

Proposition 2 The mapping ^s{a;-) defines a contraction in A^([— 1,1]) equipped with the 
Hutchinson metric. 

Clearly, ^s{(j;-) also induces a map between the Jacobi matrices associated with the related 
measures. We want now to derive the algebraic relations required to translate this mapping into 
a numerically stable procedure. 

3.2 Scaling properties of orthogonal polynomials 

Fundamental for the algebraic structure of the algorithms that we are about to develop is the 
study of the scaling properties of polynomials with respect to the maps 0. We start with the 
almost trivial, but fundamental 
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Lemma 1 LetUn{s) be any polynomial of degree n . Let {pn{ri; s)}n£'N md {pn{(T; /3)}ni£'N be the 
families of orthogonal polynomials associated with the measures 77 and a in the variables s and 
13, respectively. Then, there exists constants f^^^ , with k,r >0, such that 

Ur,m,s))^ J2 nirPk{ms)Pr{<J;P). (7) 

0<k+r<n 



Proof. The polynomial Un{4>{l3, s)) = u„((5s + (5/3) can be written in the form X]j=o n CjS-'/S" ^ . 
Expanding the monomials s^ and /3"~^ in orthogonal polynomials of and a respectively proves 
eq. 0. I 

In this paper, we shall use Lemma [1] specifically for m„(s) = pn{i]',s), the n-th orthogonal 
polynomial of the measure fj = ^s{cr; rj) in eq. (jG]). f2", with n = 0, 1, . . ., will designate uniquely 
the related coefficients and will be called the scaling matrices. For simplicity of notation, we shall 
not indicate explicitly the fj,ri,a dependence of the matrix fJ". Also, we shall let the indices k 
and r run freely, while assuming that 17)^,. = unless fc,r > and < k + r < n, so that 
effectively is a triangular matrix. The notation qn will indicate a polynomial of degree n of 
which no further specification is necessary. 

Lemma 2 For any n G N, the scaling matrix fJ" satisfies the normalization condition 

Y.^^lr? = 1. (8) 



Proof. Since / dr]{s)pn{r]] s)'^ = 1, eq. ^ implies that 

1= ((da{P)d7i{s)p^{n-A{P.s)f 



= JJda{P)d7j{s)J2 ^lr^l,r'Pk{m s)Pr{o; f3)pk' {v, s)pr' (a; 13). 

Performing the integrations with the aid of orthogonality leads to eq. ([5]). I 

Lemma 3 For any n G N, the scaling matrix extremal entries ^V^ q and i^Q „ are given by 

^n.O - -nd nZ,n = ^H^^"", (9) 

/i„(ry) ' /i„((t) 

where hn{-) denote the coefficients of the highest power in the orthonormal polynomial of degree 
n: p„{-; s) = hn{-)s" + g„-i(s). 

Proof. Clearly, because of eq. ([7]), and because Pq{(7; /3) = 1, on the one hand 

Pniv; Ss + Sf3) = rj" oP„(77; s) + q„-i(s) = 17" o/i„(?7)s" + g^_i(s), 

where <7„_i(s) and q'n-i{s) are polynomials of degree n — 1 in the variable s and the /3 dependence 
has been implied. On the other hand, 

p^ifj; Ss + 6(3) = h^{fj)S''s" + q':^,{s), 

with q'n^iis) another polynomials of degree n — 1. This proves the first eq. ([9]). Because of the 
symmetrical role of a and rj, this also proves the second. I 

Lemma 4 

s)p„(fy; 0(/3, s)) = ^ ^kASPkiv, s)Pria; /3) + ^P.((t; P)pkiv; (10) 

k,r 

where we have put Pk{-', s) := bk+i{-)pk+i{-', s) + ak( )pk{'] s) + bk{-)pk-i{-] s) and ■ can be either 
rj or a. 
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Proof. Because of eq. ([7]) 

3, s)pn{m s)) = {Ss + ~5p) ^IrPkim s)pA<y; (i) 



k,r 

and the products spkiriis) and j3pr{(j;P) can be dealt with using the recurrence relation for 
orthogonal polynomials, eq ([2]). B 

3.3 Relations between Jacobi and scaling matrices 

We can now show that the scaling properties embodied in the matrices fJ" imply algebraic 
relations between the Jacobi matrices J,j, J,^ and J^. In other words, we are capable of translating 
on these latter the action of the convolution operator: J,^ — <^j{J„] J^i), where, with a slight abuse 
of notation, we also denote by •) the induced operator in the space of pairs of Jacobi matrices. 

Lemma 5 For any n G N and S Cz [0, 1), the Jacobi matrix entry a„(?7) can be written as a linear 
combination of the Jacobi matrix entries {aj{ri), foj('7)}jLo '^"'^ ^j(''')}j=0 7 ^^th coefficients 

derived from the scaling matrix f2" ; 

an{fj) = J2 ^Ir {^[akiriWk.r + 2hk{r,)ni_^^,] + ~5[ar{a)ni^r + 2br{a)ni,,_^]) . (11) 

k.r 

The coefficients of the highest index terms an{r]) and a„(f7) at r.h.s. are (5(r25^ q)^ and (5(510 „)^, 
respectively. 

Proof. We use eq. ([2]) and orthogonality to write 

anijf) = j df]{s)spl{fj]s) (12) 
and we then expand according to the relation ([6]) and to eqs. ([7]), ((TO)) : 

aniv) = J da{(i) j dry(s)0(/3,s)[p„(77;0(/3,s))]2 = (13) 

ldcT{P)dr^{3)Y, ni^ni^^,pk>{ri; s)pr'{<J; /3)[5Pfc(ry; s)p,(a; /3) + 5P,(a; /?)pfe(r,; s)]. 

k,r,k' ,r' 

Integrations can be computed explicitly using orthogonality of the polynomials: after some ma- 
nipulations, one gets the linear relation eq. (|Tll) . Notice that the highest index with non-zero 
coefficient for both k and r is n. Direct inspection shows that the coefficient of On^rj) at r.h.s. is 
5(rj;j.o)^ and that of a„(CT) is S{^o,n?- ■ 

Lemma 6 For any n £ N, all the entries of the matrix 51""*"^ := bn+i{fj)fl"'~^^ can be computed 
linearly in terms of the entries ofil" and , with coefficients determined by the Jacobi matrix 
entries {a, {r,)}]^„ {a,{f])}]^„ {a,{a)}^^„ {bM}ho^ {bM}]=i , ^nd {b,{a)}]+^ . In fact, the 
following relation holds: 

:= br.+,{m]+' = 6 [a,{fj)n], + b,+,{^)n-^,, + b,{r^)n]_,,] + 

The entries 6„+i(ct), 6,1+1(77), only affect the computation of the extremal values 

&n+i(^)^^::ilo = sni^b^+iiv), (15) 
:= bn+iimtnii = ~snijn+i{<^). (16) 
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Proof. We start from the relation 

bn+i{fi)pn+i{ms) = {s - an{fj))pn{fj; s) - 6„ (^)p„^ i (f? ; s) , (17) 

and map s in s) = Ss + (5/3. Define Q„+i(s) := b„+i{fi)pn+i{fj] (j){l3, s)): 

Qn+i{s) = {(f){f3,s) - an{fi))p.n{fi;(j){(3,s)) - &„(?7)Pn-i('7;0(/3,s)). (18) 

Use now eqs. ([7]) and (fTO|) to express the terms at r.h.s. via the matrices fi", and Then, 
we multiply both sides of eq. ()18p by pj(r]; s)pi{a; /3) and integrate w.r.t. da{p)dri{s), to get 
eq. (dH). The r.h.s. of this equation is a linear combinations of the matrix elements of ft" and 
ft"~^. The coefficients are given by the matrix entries of J,,, Jfj and J^- Direct inspection, using 
the triangular nature of the r2"'s, reveals that the terms of highest index appearing at r.h.s. are 
0^(77), 6,1+1(77) and a„((T), &„+i((t), and it also reveals that the last part of the thesis and eqs. 
(HH) hold. I 

Lemma 7 For any n £ TSI, the Jacobi matrix entries bn+i{f]), bn+i{ri) and 6„_(-i(cr) are related 
to the scaling matrix Ct"^^ via 

k,r 

where the primed summation runs over all pairs of indices (fc, r) that are different from (n + 1, 0) 
and (0, n + 1). 

Proof. Because of eq. (|H]), 

hl+M = Y^iv^Z'f = {fi-,%\,f + {KXUf + (20) 

We can now use the explicit formulae (IT5|) and p^ . I 

4 Algorithms for the construction of the I.F.S. Jacobi Ma- 
trix 

Having derived the necessary technical lemmas in the previous section, we are now in a position 
to chain them into two algorithms for the construction of the Jacobi matrix of an affine, 
homogeneous I.F.S. described by the contraction ratio 5 and by the affine constants distribution 
a. We first develop a fixed-point, forward algorithm similar to that already exploited in [44], 
based on a technique for computing the I.F.S. convolution of two Jacobi matrices. Next, the 
solution of the fixed point equation can be obtained by a bootstrap technique like that of ref. 
|41j . This yields the closure algorithm. We will denote by J^"^ the finite truncation of any Jacobi 
matrix J to size n. 

4.1 I.F.S. convolution of Jacobi matrices 

Theorem 1 The Jacobi matrix Jf^ of the I.F.S. convolution measure f] :— ^^(cr;?]) can be com- 
puted recursively from the Jacobi matrices Ja and Jrf. Computation of a finite truncation J^"^ 
of size n of Jfj only requires knowledge of the truncated matrices J^^^ and jI^^ . 

Proof. The algorithm is structured in the following sequence of steps 
Algorithm 1. 

Computing the Jacobi matrix of an I.F.S. convolution. 

Input: the (truncated) Jacobi matrices J^""* and J^"\ the contraction factor (5, the trun- 
cation size n. 

Output: the (truncated) Jacobi matrix of fj := ^s{cr]v)- 
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0: Initialization: n = 0. One has i^Q q = 1, since po{ri;s) = po{r]:s) ~ pQ{a; l3) = 1 and 
boifj) = bo{r]) ^ bo{a) ^ 0. 

1: Induction hypothesis: {fi^, j = 0, . . . , n}, {aj{fj),j = 0, . . . , n — 1} and {bj{fj),j = 0, . . . , n} 
are known. 

2: Computation of an{ff)'- eq. in Lemma [5l 
3: Computation of the matrix $7"+^: Lemma [51 
4: Computation of bn+i{ff)'- eq. (|19p in Lemma [7] 
5: Computation of i7"+^: divide 17"+^ by 6,1+1(7^). 

6: Augment n to n + 1 and loop back to step 1 if n + 1 is less than the desired truncation n, 
otherwise stop. 

Notice that 6„+i(^) in steps [4] and [5] is never zero, if (5 > 0, or if the cardinality of the support 
of either u or ry is larger than n + 1 . B 

4.2 Fixed— point Forward Algorithm 

The previous algorithm can serve a double purpose: on the one hand, to investigate the convolution- 
like measure fy as a function of the factors a and 77. On the other hand, to obtain the Jacobi 
matrix of the invariant {5, (t)-IFS measure in an iterative fashion, as in the following: 

Algorithm 1-Fix. Computing the {&, cr) IFS Jacobi matrix. 

Input: the (truncated) Jacobi matrix ji-"'' of the distribution of affine constants cr, the 
contraction factor (5, the truncation size n and a convergence threshold. 
Output: the (truncated) Jacobi matrix J^"-* of the I.F.S. measure /i. 

0: Initialization: Let /iq be a positive probability measure in 7W([-1,1]) and let J^"^ be its 
(truncated) Jacobi matrix. 

1: For m — 1 until convergence 

Use algorithm 1 to compute Jm \ the truncated Jacobi matrix of /im := $5(17; /im-i)- 

Remark 4 Convergence of the previous algorithm is assured by the contractive nature of the 
transformation $5(0"; •), Prop. [H Numerical convergence, on the other hand, can be gauged e.g. 
according to the Frobenius norm of the difference Jm^ — J^-i- Observe that, because of theorem 
[7J this difference is exact (except for numerical errors) at any finite truncation n, that is to say, 
enlarging n does not change the values of the already computed smaller truncations of the Jacobi 
matrix. 

4.3 Closure Algorithm 

We now show that the Jacobi matrix of the invariant measure /i of a {S, cr)-IFS can be obtained 
directly and in a finite number of steps — that is, the sequence of operations of Algorithm 1-Fix 
can be exactly closed. In fact, we can set ^ = r] = fj throughout the formulae of section [3l The 
only notable difference occurs in the following lemma: 

Lemma 8 For any n, the Jacobi matrix entries bn+iifJ^) and 6„+i (cr) are related to the scaling 
matrix fl"^^^ via 

(1 - s'(-+'^)bl^M = bl^,ia)P(nu' + iZi^Z')'' (21) 

where the primed summation runs over pairs of indices (fc, r) not equal to (n-|- 1, 0) and (0, n-|- 1). 
Therefore, the Jacobi matrix entry 6„+i(/i) can be always computed from the other quantities in 
the above equation, while 6„+i(cr) is only defined from the other terms when the difference between 
the l.h.s. and the second term at r.h.s. is positive. 
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Proof. This is a consequence of Lemma [7] and can be proven similarly. I 

Remark 5 Observe the breaking of the symmetry in the role of the measures fi and a in the 
previous lemma. This is of paramount importance for the inverse problem discussed in Sect. 

Theorem 2 The Jacobi matrix of a homogeneous I.F.S. with contraction ratio S and affine 
constants distribution <j can be computed recursively from the Jacobi matrix J„ . Computation of 
a finite truncation jjf^"^ only requires knowledge of j'^P^ . 

Proof. The algorithm is effected in the following sequence of steps 
Algorithm 2. Computing the (5, cr) IFS Jacobi matrix. 

Input: the (truncated) Jacobi matrix J^"'', the contraction factor 5, the truncation size n. 
Output: the (truncated) Jacobi matrix of the I.F.S. measure /i 

0: Initialization: n = 0. One has J^q q = 1, since po(m; s) — poicr; /S) — 1 and bo{n) = &o(c) = 0- 

1: Induction hypothesis: {fl^ , j = 0, . . . ,n}, {aj(/x), j — 0, . . . , n — 1} and {6j(/x), j — 0, . . . ,n} 
are known. 

2: Computation of a„(/j,): Lemma [5] Observe that setting r] — fj — fi in eq. (ITT|) the term 
anifJ.) appears at r.h.s. with a coefficient (J(ri" q)^ ^^^^ strictly less than one. 

3: Computation of the matrix 17"+^: Lemma [51 

4: Computation of 6„+i(/i): Lemma[5] 

5: Computation of i7"+^: divide 17"+^ by 6„+i(/i). 

6: Augment n to n + 1 and loop back to step 1 if n is less than the desired truncation n, 
otherwise stop. 

Notice again that 6„+i(/i) in steps [4] and [5] is never zero, if (5 > 0, or if the cardinality of the 
support of (7 is larger than n + 1. M 

Remark 6 The algorithm can be carried out in 0{n^) operations. In addition, it can also be 
re-structured in order to compute in place of the b„{^) the squares b^di), that enter the recursion 
relation of monic orthogonal polynomials (i.e. those normalized as having unit coefficient in the 
highest power). This avoids the square root in step 4 cii^-d therefore leads to an algorithm that 
can be carried out exactly in rational arithmetics, when S and Ja are such, or symbolically, for 
instance to obtain the recursion coefficients as functions of the contraction ratio S. 

Remark 7 The storage requirement of the algorithm is 0{n^). 

Remark 8 In the classical I.F.S. case with M maps, when a is a finite sum of M atomic 
measures, eq. Q), the algorithm above can be used without any modification. In this case the 
Jacobi matrix J„ is finite and one can set bj{a) = for j > M . This entails, via eq. \14^ , that 
fc ~ for all k > M . The algorithm then runs in O^Mn?) operations and requires a storage 
of size 0{Mn), exactly as the Stieltjes algorithm in \4 _/ [ /. Indeed, the difference between the latter 
and the present algorithm is two-fold. Firstly, notice that in the polynomial pn{4>{l3j, s)), 
for any j = 1, . . . , M , is developed on the basis {pkifJ-', s)}'^^q, resulting in M coefficient vectors. 
Here, we exploit the algebraic properties of these coefficients in order to write them on the basis 
of the orthogonal polynomials of a. Secondly, the input data of the algorithms are different: in 
141^ these are the parameters of a finite set of maps — here, they are the common contraction 
coefficient d and the Jacobi matrix of the measure a. 
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4.4 Numerical examples 

Let us start by examining the fixed-point, forward algoritlim 1-Fix, choosing fiQ as the uniform, 
normahzed Lebesgue measure on [—1,1]. This choice is dictated by the simphcity of its Jacobi 
matrix, but other choices would work equahy weU. Indeed, one might even start from a suitable 
Jacobi matrix whose associated measure is not even known. On the other hand, three kind of 
possible choices for a are: 

• CTpp! a point measure with a finite number of atoms, 

• asc- a singular continuous measure, 

• aac- an absolutely continuous measure. 

We investigate these different choices because of their role in mathematical physics, as described 
in the introduction and because of the open problems on their continuity properties, discussed 
in Sect. [SJ 

Example 1 Let the atomic measure a^p ~ ^(Z?_i + Di) he made of two atoms at positions 
— 1 and 1, with equal weights. Fixing S = therefore implies that fi — Mmm^oo l^'m is the 
invariant measure of a two-maps homogeneous I.F.S. with ^ — equal weights i^i = 1^2 — \ 
and Pi = —f32 = 1. The convex hull of the support of this fractal measure is also [—1, 1]. 

Example 2 As to the second choice, asc can be taken as precisely the two-maps I.F.S. measure 
of example [71 The generated measure fi ( with a value of the contraction ratio 6 independent of 
that of Ex. Qp is now the invariant measure of an I.F.S. with uncountably many maps, whose 
fixed points are located on a Cantor set. For the case of Fig. [I] we have chosen S — 1/4. 

Example 3 Thirdly, aac can be taken as the uniform Lebesgue measure on [—1. 1]. The generated 
measure /i (with 5 is now an absolutely continuous measure supported on [—1,1], whose 

analytical properties have been studied in \50^ . where a graph of its density is also displayed. 

We exhibit the numerical convergence properties of algorithm 1-Fix, when applied to Examples 
[1] to [3l In Fig. [1] we plot, in log-linear scale, the Frobenius norms of the differences between 
the Jacobi matrices of fim and ^m-i, versus the iteration number m, computed at a fixed finite 
truncation. After an initial transient, that lasts longer for the point measure than in the other 
cases, we observe the exponential convergence typical of fixed point techniques. Obviously, the 
Jacobi matrix of the three measure just considered can be also computed via the closure algorithm 
2, with an obvious saving of computer time. In certain investigations, though, the recursive 
algorithm may be needed, when one wants to study the properties of the Jacobi matrices of a 
sequence of measures converging to the limit measure /j,. 

Let us now consider Algorithm 2. In the case when the measure a is composed of a finite 
number of atoms, both storage and cpu time can be greatly reduced, as noted above. This 
fact permits far-reaching numerical experimentations. We first consider the uniform Lebesgue 
measure on [—1,1], that can be generated by choosing 5 = 1/2 in example [1] while leaving a 
unchanged. Being the Jacobi matrix explicitly known, we have tested the error propagation, 
finding errors less that 6 • 10^^® for n as large as 250, 000, being the machine epsilon of the order 
of 2.22 10-16. 

Next, we consider the refinable functionals introduced in Ref. p8l. As a matter of facts, their 
invariant measures are supported on a different interval than [—1, 1], but our algorithms work 
equally well without change. 

Example 4 In I.F.S. language, the first numerical example discussed in Sect. 4 of Ref. 
consists of four maps, with contraction factor 5—1/2 and with fixed points l3j — j, for j = 
0, ...,3 and weights 1/8,3/8,3/8,1/8, respectively. 

We have computed the Jacobi matrix of Example H] with algorithm 2 up to n = 250,000 
without encountering any numerical instability. Because of symmetry, one can theoretically 
assess that a„ = 3/2, a value that is not automatically reproduced by the algorithm so that 
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Figure 1: Frobenius distance between the truncated Jacobi matrices of fim and ^m-i, versus m, 
for n — 4096. The three measures a described in Examples [T] to [3] are considered. The highest 
curve is for the atomic app (Ex. [U crosses), while the two lowest sets of data are for asc (Ex. O 
pluses) and aac (Ex. [3l asteriscs). 
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Figure 2: Absolute differences e„ = |a„ — 3/2| (lower data points, pluses) and e„ = |&„ — 3/4| 
(higher data points, crosses) versus n for the refinable functional defined in Example ID The 
fitting power-laws e„ ~ ri^ have exponent 7 — .85 and 7 = —2, respectively 

it can be used to gauge numerical error propagation. In figure [2] we observe a mild, slower 
than linear error growth for the diagonal coefficients, that is orders of magnitude lower than 
the difference between 6„ and the asymptotic limit boo — 3/4, also reported in the figure. The 
observed power-law decay of such difference, with exponent 7 = —2, is therefore to be deemed 
reliable and suggestive of the absolute continuity of the orthogonality measure fj,, according to 
the theory presented in Sect. [5j 

To provide a more stringent verification, we have also computed, in the same case as above, 
the difference between the results of Algorithm 2 and those of the Stieltjes technique of Ref. [H], 
this latter run in quadruple precision (machine epsilon of the order of 10"'^^). The slower than 
linear growth of the absolute differences in the diagonal and the outdiagonal components of the 
Jacobi matrix, reported in Fig. |3l suggests the numerical stability of both techniques. 

Finally, we have run the same test on the second example proposed by Laurie [38], who 
observed a linear growth of the error in his technique, discovering a comparable behavior, on the 
larger range here displayed: Ex. [5] and Fig. |4l Using the theory of the next section, absolute 
continuity of the associated measure can also be conjectured from the computed J^. 
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Figure 3: Absolute differences e„ = la„ — a„| (crosses) and e„ = |6„ — b\ (pluses) versus n between 
the a„, bn provided by algorithm 2 and the a„, 6„ from the algorithm of Ref. [41] when applied 
to the same case of Ex. IH Fig. [2] 
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Figure 4: Absolute differences e„ ~ |a„ — a„| (crosses) and e„ = |6„ — 6| (pluses) versus n between 
the a„, bn provided by algorithm 2 and the a„, 6„ from the algorithm of Ref. |41j when applied 

to Ex. m 

Example 5 In example^^ use the modified weights 1/8,3/8,3/8,1/8. 

We propose numerical experiments on the case of infinite Jacobi matrices Ja later on, in Sect. 
[71 We now pause for a theoretical digression of some interest. 

5 Analytical properties of the invariant measure 

The theory developed in the previous section can be also employed for an ambitious goal: to 
study numerically the analytical properties of the invariant measure /i of a {S, cr)-IFS. To give an 
idea of what we believe can be achieved along this line, in this section we briefly discuss two types 
of results, one for measures whose support is a full interval, the other for measured supported on 
Cantor sets. 

The continuity properties of the measure fi follow in a complicated way from those of a and 
from the contraction ratio S. For instance, even in the realm of conventional, finite cardinality 
I.F.S., cases of point measures a leading to either singular continuous, or absolutely continuous 
measures are well known. Reference [SO] is an attempt to attack this problem in full generality: 
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Figure 5: Differences |6n — ^| (decreasing curves) and partial sums Sn (increasing curves) versus 
n for the three measures described in the text. Inside each group the curves arrange themselves 
for large n from bottom to top, starting from Ex. [I] with S — Si, followed hy 5 = 82 and finally 
S = 53. 



typically, we find that fi is "more continuous" than a. We have proven that when a is absolutely 
continuous with a bounded density, so is fi, for any d. This covers most cases commonly encoun- 
tered in numerical analysis, but not the general situation that we discuss in this paper. Various 
techniques have been proposed in [50] to verify numerically the continuity type of a measure /x. 
We now present a different one. 

Recall that the Nevai class of measures N{aoo, boo) contains the orthogonality measures asso- 
ciated with Jacobi matrices for which a„ Ooo and 6„ — >■ 600, as n — ^ 00. From the speed of this 
convergence one can infer the continuity properties of /i (see e.g. |61].|10]). For instance, when 

^ |an - acx)| + |6n - feool < 00 (22) 

n 

the measure fi is absolutely continuous w.r.t. the Lebesgue measure on the interval [aoo — 

2^00; ^oo ~t~ 2600]. 

The numerical stability and the performance featured by Algorithm 2 permit to compute 
Jacobi matrices of large orders and therefore to hint to the existence of Ooo and boo and to the 
validity of eq. (|22|) . We now apply this technique to the Erdos problem of infinite Bernoulli 
convolutions. For details of this problem see the review paper |53j and references therein. 

Example 6 Let a — ^{Dq + Di) (a minor variation of example\T\) and select three values of 5, 
61 = 2-1/2 _ 0.7071067811865, 62 = 3/4 = .75 and 63 = 1/pi - 0.7548776662467, pi being a 
Pisot number. 

It is rigorously known that the generated measure ^ is absolutely continuous in the first case 
and singular continuous in the third. It is also absolutely continuous for Lebesgue almost all 
values of 5 between one half and one, but it is only conjectured that rational values in this 
interval, such as 82, belong to this case. Notice that 82 is very close to ^3. Also notice that 
''oo = 1/2 is exactly known in the three cases. 

In figure[5]we plot both |6„ — ^| and the partial sums Sn J2j ~ 1/^1 versus n in doubly 
logarithmic scale for the three cases listed. Convergence of Sn is observed in the first two cases, 
as can be inferred by the exponent of the power-law decay of |5„ — ^|. In the third case, this 
decay (which is also evident) is too slow to imply the absolute continuity of fi. Therefore, our 
technique gives results that are consistent with rigorous facts and, which is more important, with 
the conjectured absolute continuity for the case S = 62- 

The measures just discussed are supported on the full interval [0, 1]. Suppose now that 5*^ is 
a Cantor set. The measure fi is then singular continuous and does not belong to a Nevai class. 
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Our theory permits to derive estimates for the capacity of S*^, in the sense of potential theory 
|59| . Precisely, we have 

Proposition 3 Suppose that the distribution of fixed points a and the invariant measure /i are 
regular, in the sense of potential theory. Let then and be the capacities of the supports of 
a and of ^i, respectively. Then, 

C,<C^<C„+ log(ri) <C, + 5. (23) 



Proof. Consider the coefhcient hn{ii) of s" inp„(/i;s): from eq. ^ it follows that 

n 

h^{^Ji) = l/X{h,{^i). (24) 

If jji is any regular probability measure, the asymptotic relation hn[n) ~ 6"*-^^ holds [59]. Then, 
l^"!^! behaves asymptotically as e"^, where A = — Co-. Since the support of a is enclosed in 
that of 11 this difference is always positive. Since < 1 for all k,r, using Lemma [H we get 

/l„(^) < ft.„(cr)^~", 

and the result follows. I 

Remark 9 It is possible to prove the regularity of I.F.S. singular continuous measures fi in large 
generality, see J59f . 

Observe finally that computing large order Jacobi matrices can lead to a numerical estimate 
of the capacity of the support of /Lt, via the quantity i log(ft.„(/i)) = —^Y^^=i^og{bj{fi)) that 
converges to C^. Convergence is typically slow, but it can be accelerated by suitable extrapolation 
techniques [8]. 



6 A spectral data technique 

The numerical determination of the Jacobi matrix of fj in Sect. [U eq. ^ can also be effected 
without recurring to the algebraic theory developed in Sect. [3] and leading to algorithms 1 and 
2. In fact, the problem can be cast into a forward/inverse Gaussian quadrature determination. 
Nonetheless, this approach only provides us with an analogue of the recursive algorithm 1-Fix 
and not of the faster algorithm 2. The key observation is the following: 

Lemma 9 The formula 

n n 

rf^(s)/(5) = ^5]^"^(r7)u;i")(a)/(<s4")(r;) + 54")(a)), (25) 
j=i k=i 

where x^"'*(-) and Wj^\-) are Gaussian points and weights, respectively, is exact for f G P2n-i- 
Therefore, it can be used to compute J^"'' exactly. 

Proof. Let / in eq. ([6]) be a polynomial of degree at most 2n — 1 in the variable s. Then, 
f{4){j3,s)) — f{5s + SP) can be exactly integrated with respect to the measure r] by n-points 
Gaussian summation in the variable s. This Gaussian formula can be easily obtained by the 
spectral problem of Also, f{4>{l3,s)) is a polynomial of degree 2ri — 1 in the variable /3, 

and the same is its integral w.r.t. dri{s). This latter polynomial can be exactly integrated with 
respect to the measure a by an n-points Gaussian summation obtained from J^""* . I 

The above Lemma can also serve as an alternative proof of Thm. [TJ Next, observe that the 
r.h.s. of eq. (|25|) is the integral of / with respect to the sum of atomic measures. Stable 
algorithms for computing the Jacobi matrix of a finite sum of atomic measures, due among others 
to De Boor and Golub [TT], Gragg and Harrod [21], Fischer [50], Reichel [53] and Laurie [37] are 
well known and can be put to use, to give the following: 
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Algorithm 3. Computing the I.F.S. convolution. 

Input: the (truncated) Jacobi matrices ji-"'' and the contraction factor 5, the trun- 
cation size n. 

Output: the (truncated) Jacobi matrix of fj $5(cr;?7). 
1: Compute Gaussian points and weights for a from J^-"'. 

(n) 

2: Compute Gaussian points and weights for r/ from . 

3: Using Lemma [9] compute J^"'' using one of the algorithms just quoted. 

It is immediate to obtain an iterative version, 3-Fix, along the same lines of Alg. 1-Fix: 
from step [3] loop back to step [2] replacing rj by fj. This fixed-point algorithm then provides us 
with the Jacobi matrix of ^. As a matter of facts, algorithm 3-Fix works fine as far as absolutely 
continuous measures a (like that in Example[2]) are involved. Instead, when /i {not a) is supported 
on a Cantor set (like e.g. in Example [1]), we have observed that it achieves convergence only 
for Jacobi matrices of the size of about a thousand. This is due to the fact that the n? points 
in Lemma IHl crowd around a fractal and the relative precision in their distance diminishes. As 
explained in detail by Laurie [211 (see also [211 HO] ) this fact impairs the reconstruction of the 
Jacobi matrix from the Gaussian points and weights. One can therefore appreciate by comparison 
the computational advantage brought about by the algebraic theory of Sect. [3] 

7 I.F.S. Quadratures and the Inverse Problem 

The algorithm 2 presented above can be reversed, in order to compute from J^. This is the 
basis of the solution of an inverse problem, that can be used in an approximation problem: that 
of finding I.F.S. quadratures [71 149) . 

Definition 3 Given a target measure fi, whose support is enclosed in a finite interval, an I.F.S. 
quadrature for fi is a sequence of I.F.S. measures fi''"^ that satisfies — for any n £ N. 

Remark 10 Def. ^ implies that /i^"^ integrates exactly polynomials up to degree 2n — 1 and 
therefore the sequence {/i^"^} is weakly convergent to /i. Clearly, the linear combination, with 
Gaussian weights, of the atomic measures sitting at the Gaussian points of order n is an I.F.S. 
quadrature, degenerate in the sense that Gaussian points are the fixed points of a finite set of 
maps with contraction rate 6 = 0. 

Def. @ formalizes a truncated inverse problem, in the family of fractal inverse problems 
[1 [3 [71 ISHU]. If we now restrict ourselves to (5, (t)-I.F.S. of the kind ®, the following 
approximation result guarantees that solutions do exist: 

Theorem 3 ([33]) Let fi he a measure with an infinite number of points of increase. Then, for 
all n > there exists SnifJ-) > such that for all 6 G [0,(5„(/i)) there exists an homogeneous affine 
I.F.S. with n maps that satisfies Def. 

This theorem and Def. [Slalso imply that any finite symmetric tridiagonal matrix with positive 
out-diagonals is the truncation of the Jacobi matrix of a (^, ct)-I.F.S. with non-vanishing S. 
We shall now apply this theorem while developing the inverse of algorithm 2, in a form that is 
also suitable for the numerical determination of the maximal value ^„(/x). 

Theorem 4 ([49]) The truncated Jacobi matrix of the distribution of fixed points a of a 
{6,a)-I.F.S. with contraction ratio d that provides an I.F.S. quadrature of a measure fi can be 
computed recursively from the truncated Jacobi matrix JJ^, provided 6 < 5„(/x). Whether the last 
condition holds can be verified recursively. 

Proof. The algorithm is effected in the following sequence of steps 
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Algorithm 4. Solving the inverse I.F.S. problem. 

Input: the (truncated) Jacobi matrix J^"-* of the target measure the contraction ratio 
i5, the maximum size n. 

Output: the (truncated) Jacobi matrix of u, the largest allowed truncation size n . 

0: Initialization: n = 0. One has ilg o — li since Po(m; s) = pQ{a; /3) = 1, and ba{p) = bo{a) = 
0. 

1: Induction hypothesis: {fi-',j = 0, . . . ,n}, {aj{a),j = 0, . . . ,n — l}, and {bj{a),j = 0, . . . ,n} 
are known. 

2: Computation of a„(a): Lemma [5] Observe that a„((7) has a non-zero coefficient in eq. 
pT|) . due to Lemma 131 

3: Computation of the matrix $7"+^: Lemma [51 

4: Computation of 6^^]^ ((t): Lemma [51 eq. (|^T|) . 

4: Stopping alternative: either &^_|_i(cr) > 0, therefore continue, or else 6 > 5„+i(/i), h — n 
and stop. 

6: Computation of £7"+^: divide 17"+^ by 6„+i(/i). 

7: If n < n augment n to n + 1 and loop back to 1, else h — n and stop. 

Remark 11 When termination occurs at a certain value of n = n < h at step 4-, then 5 is larger 
than Sn+i{fi), but smaller than (5„(/i). Therefore, using Algorithm 4 iteratively at different values 
of 5, one can determine the sequence of values 6n, at varying n. In principle, the algorithm never 
stops only if the target measure fj, is exactly generated by an affine IFS with contraction ratio S 
and a measure a with an infinite number of points of increase. 

To establish the numerical stability of both the forward algorithm 2 and the reverse algorithm 
4 we have chosen a target Jacobi matrix of particular significance, the Fibonacci tridiagonal 
matrix, whose orthogonality measure is singular continuous. This measure is not the invariant 
measure of a (5-homogeneous I.F.S. and yet, as seen above, any finite truncation of its Jacobi 
matrix coincides with the truncation of the Jacobi matrix of a ((5, (t)-I.F.S. 

Example 7 Let a„ = for all n, and let bn take either the value A = 2/5 or the value B = 1/2. 
These values are arranged in the aperiodic Fibonacci sequence ABABBA..., generated by the 
substitution rules A — > AB, B A on the seed A. It has the property that two A 's never follow 
each other, but are separated by at most two B's. 

Using algorithm 4 we have computed the sequence (5„(/.t) at increasing values of n up to n, 
as well as the Jacobi matrix J^""* for a feasible value of 5, close to the maximum allowed value 
6n{ii). We have then applied Algorithm 2 to recompute the original Fibonacci Jacobi matrix. 
While the null diagonal entries are recovered exactly because of the nature of the algorithms. 
Figure [6l plots the absolute errors e„ in the reconstruction of the sequence of 6„. The observed 
behavior, confirmed by other experiments, is the most convincing experimental verification of 
the numerical stability of the direct and inverse algorithms 2 and 4 presented in this paper. 

8 Conclusions 

We have presented a new family of algorithms for the direct/inverse computation of the Jacobi 
matrix of the invariant measure of a homogeneous afSne I.F.S.. Experimental results suggest 
that these algorithms are stable to large orders when programmed in floating point arithmetics. 

On the one hand, this remarkable stability calls for a detailed error analysis, that should 
unveil the reasons why the algebraic treatment of Sect. [3lis more stable than any other existing 
technique (to the author knowledge) and in particular than the Gaussian technique of Sect. [6l 
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Figure 6: Example [T] Absolute errors e„ versus n, for the reconstruction of the out-diagonal 
Fibonacci Jacobi matrix coefficients 6„, a.t 6 — 1.119837 10^^, the maximum allowed value for 6 
at n = 3500 being computed as approximately (535oo(m) = 1-124611 10^^. 

We conjecture that this is due to the fact that in our approach we exclusively deal with Jacobi 
matrices, but a more thorough investigation, that we plan to develop in further publications, is 
in order. 

On the other hand, the versatile tools that have been introduced in this paper can now be 
applied to a variety of problems, both from the theoretical and from the applied side. In the first 
respect, we would like to investigate to what extent we can infer the fine structure properties 
of the generated measure fj, from those of a and whether more can be said from the potential 
theoretical point of view, e.g. on the asymptotic properties of the sequence of Jacobi matrix 
entries and on the Fourier transform of /i and of its orthogonal polynomials. In the second 
respect, we would like to evaluate the full potential of the approximation/inverse problem of 
Sect. [7] on significant problems, until now beyond the reach of conventional algorithms. 

Acknowledgements It is a pleasure to have this opportunity to thank Dirk Laurie for providing 
his code for [37j and for related interesting discussions. 

References 

[1] S. Abenda, G. Turchetti, Inverse problem for fractal sets on the real line via the moment 
method, Nuovo Cim. B 104, 213-227 (1989). 

[2] N.I. Akhiezer, The Classical Moment Problem, Hafner, New York, NY. (1965). 

[3] M.F. Barnsley and S.G. Demko, Iterated function systems and the global construction of 
fractals, Proc. R. Soc. London A 399, 243-275 (1985). 

[4] M. F. Barnsley, V. Ervin, D. Hardin, and J. Lancaster, Solution of an inverse problem for 
fractals and other sets, Proc. Natl. Acad. Sci. U.S.A. 83, 1975-1977 (1986). 

[5] M. F. Barnsley, Fractals Everywhere, Academic Press, New York, NY (1988). 

[6] B. Beckermann, E. Bourreau, How to choose modified moments?, J. Comput. Appl. Math. 
98, 81-98 (1998). 

[7] D. Bessis, S. Demko, Stable recovery of fractal measures by polynomial sampling, Physica D 
47, 429-438 (1991). 

[8] C. Brezinski, M. Redivo Zaglia, Extrapolation Methods: Theory and Practice. North Holland, 
Amsterdam (1991). 



17 



[9] C.W. Clenshaw, A Note on the Summation of Chebyshev Series, Math. Tables Aids Comput. 
9, 118-120 (1955). 

[10] D. Damanik and B. Simon, Jost function and Jost solution for Jacobi matrices, I, Invent. 
Math. 165, 1-50 (2006). 

[11] C. de Boor and G.H. Golub, The numerically stable reconstruction of a Jacobi matrix from 
spectral data. Linear Alg. Appl. 21, 245-260 (1978). 

[12] S. G. Demko, Euler Maclauren type expansions for some fractal measures, in H.O. Peitgen, 
J.M. Henriques and L.F. Penedo Eds. Fractals in the Fundamental and Applied Sciences, 
Elsevier-North Holland, Amsterdam pp. 101-110 (1991). 

[13] P. Diaconis, M. Shahshahani, Products of Random Matrices and Computer image Genera- 
tion, Contemporary Math., 50, 173-182 (1986). 

[14] E. Diekema, T.H. Koornwinder, Differentiation by integration using orthogonal polynomials, 
a survey. larXiv:1102. 52191 ^1 [math.CA], (2011). 

[15] G. Donovan, J. Geronimo, D. Hardin and P. Massopust, Construction of orthogonal wavelets 
using fractal interpolation functions, SIAM J. Math. Anal. 27, 1158-1192 (1996). 

[16] J. H. Elton and Z. Yan, Approximation of measures by Markov processes and homogeneous 
affine iterated function systems. Constr. Appr., 5, 69-87 (1989). 

[17] C. Escribano, A. Giraldo, M.A. Sastre, E. Torrano, Computing the Hessenberg matrix as- 
sociated with a self-similar measure, J. App. Theory 163, 49-64 (2011). 

[18] H.-J. Fischer, On the Condition of Orthogonal Polynomials via Modified Moments , Z. Anal. 
Anwendungen 15, 223-244 (1996). 

[19] H.-J. Fischer, Recurrence Coefficients of Orthogonal Polynomials with Respect to Some 
Self-Similar Singular Distributions, Z. Anal. Anwendungen 14, 141-155 (1995). 

[20] H.-J. Fischer, On generating orthogonal polynomials for discrete measures, Z. Anal. Anwen- 
dungen 17, 183- 205 (1998). 

[21] B. Forte, E.R. Vrscay, Solving the Inverse Problem for Measures Using Iterated Function 
Systems: A New Approach, Adv. Appl. Prob. 27, 800-820 (1995). 

[22] W. Gautschi, On the Construction of Gaussian Quadrature Rules from Modified Moments, 
Math. Comp. 24, 245-260 (1970). 

[23] W. Gautschi, Computational Aspects of Orthogonal Polynomials, in Orthogonal Polyno- 
mials, P. Nevai Ed., Kluwer, Dordrecht pp. 181-216 (1990). 

[24] W. Gautschi, Orthogonal polynomials: computation and approximation. Numerical Mathe- 
matics and Scientific Computation, Oxford University Press, New York (2004). 

[25] W. Gautschi, Orthogonal polynomials (inMatlab). J. Comput. Appl. Math. 178, 215-234 
(2005). 

[26] W. Gautschi, On generating orthogonal polynomials, SIAM J. Sci. Comp. 3, 289-317 (1982). 

[27] W. Gautschi, L. Gori, F. Pitolli, Gauss quadrature for refinable weight functions, Appl. 
Comp. Harm. Anal. 8, 249-257 (2000). 

[28] G.H Golub and J.H. Welsch, Calculation of Gauss Quadrature Rules, Math. Comp. 23, 
221-230 (1969). 

[29] G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature rules. Math. Comp., 23, 
221-230 (1969). 



18 



[30] G. H. Golub and G. Meurant, Matrices, moments and quadrature with applications, Prince- 
ton Univ. Press (2010). 

[31] W. B. Gragg and W. J. Harrod, The numerically stable reconstruction of Jacobi matrices 
from spectral data Numer. Math. 44, 317-335 (1984). 

[32] I. Guarncri and G. Mantica, Multifractal Energy Spectra and their Dynamical Implications, 
Phys. Rev. Lett. 73 3379-3382 (1994). 

[33] C.R. Handy and G. Mantica, Inverse Problems in Fractal Construction: Moment Method 
Solution, Physica D 43 17 36 (1990). 

[34] J. Hutchinson, Fractals and self-similarity, Indiana J. Math. 30, 713-747 (1981). 

[35] P. Janardhan, D. Rosenblum and R. S. Strichartz, Numerical experiments in Fourier asymp- 
totics of Cantor measures and wavelets. Experiment. Math. 1, 249 273 (1992). 

[36] D.P. Laurie, Computation of Gauss-type quadrature formulas. J. Comput. Appl. Math. 127, 
201-217 (2001). 

[37] D. Laurie, Accurate recovery of recursion coefficients from Gaussian quadrature formulae, 
J. Comp. Appl. Math. 112, 165-180 (1999). 

[38] D. Laurie and J. De Villiers, Orthogonal polynomials and Gaussian quadrature for refinable 
weight functions, Appl. Comp. Harm. Anal. (2004) 17, 241-248 (2004). 

[39] D. Laurie and J. De Villiers, Orthogonal polynomials for refinable linear functionals. Math. 
Comp. 75, 1891 1903 (2006). 

[40] D.P. O'Leary, Z. Strakos, P. Tichy On the sensitivity of Gauss-Christoffel quadrature, Nu- 
mer. Math. 107, 147-174 (2007). 

[41] G. Mantica, A Stieltjes Technique for Computing Jacobi Matrices Associated With Singular 
Measures, Constr. Appr.,12, 509-530 (1996). 

[42] G. Mantica, Quantum intermittency in almost periodic systems derived from their spectral 
properties, Physica D, 103, 576-589 (1997). 

[43] G. Mantica, Wave propagation in almost-periodic structures, Physica D, 109, 113-127 

(1997). 

[44] G. Mantica, On Computing Jacobi Matrices associated with Recurrent and Mobius Iterated 
Functions Systems, J. Comp. Appl. Math., 115, 419-431 (2000). 

[45] G. Mantica, Fourier Transforms of Orthogonal Polynomials of Singular Continuous Spectral 
Measures, Int. Ser. Numer. Mathematics 131, 153-163 (1999). 

[46] G. Mantica, S. Vaienti, The asymptotic behaviour of the Fourier transform of orthogonal 
polynomials I: Mellin transform techniques, Ann. Henri Poincare 8, 265-300 (2007). 

[47] G. Mantica, D. Guzzetti, The asymptotic behaviour of the Fourier transform of orthogonal 
polynomials II: Iterated Function Systems and Quantum Mechanics, Ann. Henri Poincare 8, 
301-336 (2007). 

[48] G. Mantica, Fourier-Bessel functions of singular continuous measures and their many asymp- 
totics, Electron. Trans. Numer. Anal. (Electronic), 25, 409-430 (2006). 

[49] G. Mantica, Polynomial Sampling and Fractal Measures: I.F.S.-Gaussian Integration, Num. 
Alg. (2007) 45 269-281. 

[50] G. Mantica, Dynamical Systems and Numerical Analysis: the Study of Measures generated 
by Uncountable I.F.S, Num. Alg. (2010) 55 321-335. 



19 



[51] D. Mauldin and M. Urbansky, Dimensions and measures in infinite iterated function systems, 
Proc. London Math. Soc. (3) 73, 105-154 (1996). 

[52] F. Mendivil, A generalization of IFS with probabilities to infinitely many maps. Rocky 
Mountain J. Math. 28, 1043-1051 (1998). 

[53] Y. Peres, W. Schlag, B. Solomyak, Sixty Years of Bc;rnoulli Convolutions, in Fractal Geom- 
etry and Stochastics II, Progress in Probability 46 39-68, (2000) Birkhauser, Basel. 

[54] L. Reichel, Construction of polynomials that are orthogonal with respect to a discrete bilinear 
form, Adv. Comp. Math. 1, 241-258 (1993). 

[55] R.A. Sack and A.F. Donovan, An Algorithm for Gaussian Quadrature Given Generalized 
Moments, Dept. of Maths Publ., Univ. of Salford, England (1969). 

[56] R. Strichartz, Self-similar measures and their Fourier transforms 1, Indiana U. Math. J. 39, 
797-817 (1990). 

[57] R. Strichartz, Self-similar measures and their Fourier transforms II, Trans. Amer. Math. 
Soc. 336, 335-361 (1993). 

[58] R. Strichartz, Self similar measures and their Fourier transforms III, Indiana University 
Mathematics Journal, 42 367-411 (1993). 

[59] H. Stahl, V. Totik, General Orthogonal Polynomials, Cambridge University Press, Cam- 
bridge (2010). 

[60] W. Sweldens and R. Piessens, Quadrature formulae and asymptotic error estimates for 
wavelet approximation of smooth functions, SIAM J. Numer. Anal. 31, 1240-1264 (1994). 

[61] W. Van Assche, Asymptotics for orthogonal polynomials and three-term recurrences, in 
Orthogonal Polynomials; Theory and Practice NATO-ASI series C issue 294, pp. 435-462 
(1990). 



20 



